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SUMMARY 


A generalized coordinates form of the penalty finite element algorithm 
for the three-dimensional parabolic Navier-Stokes equations for turbulent 
subsonic flows has been derived. This algorithm formulation requires only 
three distinct hypermatrices, and is applicable using any boundary fitted 
coordinate transformation procedure. The tensor matrix product approxima- 
tion to the Jacobian of the Newton linear algebra matrix statement has 
been derived. The Newton algorithm has been restructured to replace large 
sparse matrix solution procedures with grid sweeping, using ex-block 
tridiagonal matrices, where a equals the number of dependent variables. 
Numerical experiments have been conducted, and the resultant data gives 
guidance on potentially preferred tensor product constructions for the 
penalty finite element 3DPNS algorithm. 
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INTRODUCTION 


The time-averaged, three dimensional Navier-Stokes equations, for 
steady, turbulent, subsonic flow of a compressible, heat conducting fluid, 
can be simplified to admit use of an efficient space marching numerical 
solution procedure under certain restrictions. Baker, et.al. [1] documents 
the derivation of the simplified, so-called "parabolic Navier-Stokes" 
equation set using formal ordering arguments. References [2-7] document 
application of a penalty, finite numerical solution algorithm for the 
parabolic Navier-Stokes equations for a variety of subsonic flow configu- 
rations, including an embedding within an interaction algorithm to impose 
axial pressure gradient feedback. These results have provided a basic 
assessment of the accuracy, convergence and versatility aspects of the 
finite element penalty algorithm, as well as detailed comparison between 
experimental data and prediction for various turbulent flow geometries. 

These numerical predictions have been accomplished using the CMC:3DPNS 
computer code [8-10] which has evolved and matured principally under 
support of this contract. With the theoretical construction well verified, 
the requirement to improve efficiency becomes of next importance. The 
specific goals of this contract modification were to construct and evaluate 
a matrix tensor product factorization of the Newton iteration algorithm and 
to investigate inclusion of an embedded Poisson equation solution procedure. 
The results of pertinent analyses are documented in this report. 
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SYMBOLS 


a 

A 

b 

B 

C 

Eu 

f 

F 

h 

H 

i 

j 

J 

k 

«. 

L 

M 

n 

N 

P 

q 

Q 

R n 

Re 

s,S 


U 

x i 

8 

9R 

6 

6Q 


boundary condition coefficient; parameter 
initial value matrix; hypermatrix prefix 
constant 

hypermatrix prefix; matrix 
turbulence model coefficient 
Euler Number 

function of known argument 

finite element matrix; discretized equation system 

metric coefficient 

stagnation enthalpy; Hilbert space 

i ndex 

i ndex 

Jacobian 

finite element basis degree; turbulence kinetic energy 

summation index; differential operator 

differential operator 

number of finite elements spanning R n 

unit normal vector; dimension of space 

finite element cardinal basis; discrete index 

pressure; iteration index 

heat flux vector; generalized dependent variable 
generalized semi -discrete dependent variable 
spatial domain of differential operator 
Reynolds Number 
source term 

finite element assembly operator 
velocity vector 

Reynolds kinematic stress tensor 
convection matrix 
Cartesian coordinate system 
partial derivative operator 
boundary of solution domain R n 
Kronecker delta; parameter 
iteration vector 
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A mesh measure; increment 

c isotropic dissipation function; parameter 

C transformed coordinate 

n.j curvilinear coordinate system 

k heat conductivity coefficient 

A multiplier 

v kinematic viscosity 

p density 

Qjj Stokes stress tensor 

£ summation 

<}> constraint dependent variable 

w sublayer damping function 

ft solution domain 


Superscripts 


e 

h 

P 

T 


finite element reference 
solution approximation 
iteration index 
matrix transpose 
ordinary derivative 


Subscripts 


a 

e 

i > j > k , £ 

j 

o 


dependent variable index 
finite element reference 
tensor indices 
integration step index 
reference state 


Notation 

{ ) 

[ ] 
u 


column matrix 
square matrix 
union 

belongs to 
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PROBLEM STATEMENT 


Overview 

The basic requirements of a numerical algorithm construction for the 

three-dimensional parabolic Navier-Stokes (3DPNS) equation set are accuracy, 

geometric versatility and efficiency. Any theoretical basis, e.g., finite 

difference, finite volume, finite element, etc., applied to constructing a 

3DPNS algorithm ultimately yields the linear algebra statement {F} = {0}, 

where elements of { F } are strongly nonlinear functions of the dependent 

variable set q (x.). Dependent upon the algorithm designer's decisions, 
cc j 

members of the set q can include density, mean velocity vector, stagnation 

ot> 

enthalapy, turbulent kinetic energy, isotropic dissipation function, pressure 
and/or scalar potential fields and six components of the Reynolds stress 
tensor. 

For the subject finite element penalty algorithm statement, and the 
3DPNS equation set, the functional form of the algebra statement is, 

(FI (k, A,e, AX!, {QI }) } = {0} (1) 

In equation 1, k is the polynomial degree of the finite element basis 

h 

selected to construct the semi-discrete approximation q n ( x - ) to q (x.), 

ot J CX J 

X is the penalty function scalar multiplier, 0 is the implicitness factor 
of the downstream integration step Axj, where Q = H is the trapezoidal 
rule, and {QI ( x i ) } is the array of expansion coefficients of q^(xj), evaluated 
at the node coordinates of the (spatial) discretization UR 2 of the 3DPNS 
solution domain ft = R 2 x Xi. For all 0 > 0, equation 1 is a nonlinear 
algebraic equation system eligible for solution using any of a multitude of 
(approximation) procedures. These all can be interpreted within the frame- 
work of the basic Newton iteration algorithm matrix solution statement. 

cj( IQ 1 >)] j+i («on!$ - - (2) 

where p is the iteration index at step x j + -| * ar, d 

«»$t! * + Wljt] 


(3) 



The Jacobian [J] appearing in equation 2, is by definition the square matrix, 

; |m> (4 

where both I and J range 1 16. 

This contractual project phase initiated evaluation of decisions made in 
construction of equations 1-4, as they principally affect algorithm 
accuracy and efficiency. 

Parabolic Navier-Stokes Equations 

The three-dimensional parabolic Navier-Stokes (3DPNS) equations are a 
simplification of the steady, three-dimensional time-averaged Navier-Stokes 
equations, which in Cartesian tensor conservation form are 


L(p) = 


9 Xj 


pu, 


(5) 


L K> - lx: 
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i u j + p6 ij + p Vj - 


°ij 
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In equations 5-9 the usual superscript bar notation denoting time-averaged 
quantities [11] has been deleted for clarity. The time-averaged dependent 
variables are density (p), mean momentum vector (pu^ ), pressure (p) and 
stagnation enthalpy (H). Further, 5.. denotes the Kronecker delta, and the 

* J 

Stoke's stress tensor (o..) and heat flux vector (q^ are defined as, 

* \J 





( 10 ) 



where p and k are laminar viscosity and heat conductivity respectively. 
E-. is the symmetric mean flow strain rate tensor. 


(ID 





( 12 ) 


Finally, -u^u . is the symmetric Reynolds stress tensor with trace equal 

to 2k, where k is the turbulent kinetic energy. For present purposes, 

uTuT is assumed correlated in terms of k, u. and e, the isotropic dissipa- 
* J J 

tion function, in the form. 


Vj 


e c.ks. . 

i ij 



C ’ C ^ 1k E kJ 


(13) 


where the C a , 1 < a _< 4, are known constants [1]. 

The parabolic approximation to the steady flow Navier-Stokes set, equations 
5-13, is generated by assuming a principal direction of the flow persists, 
say parallel to the (curvi-1 inear) coordinate xi. Assuming the corresponding 
mean velocity component ui is of order unity, i.e., 0(1), and that the other 
two orthogonal components u^ are smaller, say of 0(6), then for modest 
density variation the continuity equation 1 confirms that for 3/3xi = 0(1), 
then 3/3x^ = 0(6" 1 ). Proceeding through the analysis details [12] confirms 
that the 0(1) 3DPNS equation set is 
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(14) 



= 0 
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pU5>U £ ‘ 0l £ 


= 0 


05 ) 


L ^ pl V 3x. 
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p5 k« * pu k u i 


= 0 


( 16 ) 




pHu . 
J 


+ T ~ 
9x. 


pH'u 


3H 

£ K 3x, 


Vi£ - U i°i£ 


= 0 (17) 


L(pk) = 


3x 


pku . 


3x, 


< P ‘kT U ?£- »*<o) 9k 


jr 3x . 
J 


3u ] 


+ w + pc a 0 
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(18) 


L(pc) = 


_ _3_ 

3x . 
J 


Pk,U. 


3x, 


pc -4 u:ur 3e 


e e j £ 3x 


J 


3ui 

+ c e pu ru Ti 


+ C ^pc 2 /k=0 ( 19) 

In equations 14-19 the tensor index summation convention is 1 < (i,j) < 3 
and 2 < (k, £) <_ 3. The Reynolds stress tensor constitutive equation 13 
also becomes considerably simplified under the ordering analysis. For example, 
in rectangular Cartesian coordinates, and retaining the first two orders of 

terms yields equation 20. An equation of state p = p(p,H) completes the basic 
3DPNS statement. 


\ 

\ 
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3DPNS Equation Set Completion 

Equations 14-20 do not represent a well -posed, initial -boundary value 
problem statement for the dependent variable set. In particular, equation 
14 is the sole definition for u^ , yielding an underdetermined system. The 
finite element penalty algorithm construction of Baker [1] yields a well- 
posed problem statement by inclusion of both 0(6) transverse momentum 
equations, definition of an auxiliary harmonic function for a penalty 
constraint, and definition of complementary and particular solutions to a 
pressure Poisson equation formed from equation 16. 


The derived pressure Poisson equation is, 


L (p) s- 


3L(pu j? ) 


3x, 


= Jfj> 

K 


3 V*k 


IpVi 


L 


0 


( 21 ) 



the solution to which is defined as 


P( x f ) = P c (x £ , x,) + p p (x £ , x,) 


( 22 ) 


The complementary pressure p c (x £ , x,) is the solution to the homogeneous 
form of equation 21 with exterior flow boundary conditions [1], The 
particular pressure is computed throughout the 3DPNS solution from equation 
21, and added in a delayed manner into equation 22, used in the u, momentum 
equation solution, yielding a multi-pass interaction algorithm. 

The retained 0(d) transverse momentum equation set is 


3 

- 

3 

r i 

3x . 
J 

P Vj_ 

3x k 

Q 

7T 

Xd 
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The auxiliary harmonic variable <f> is defined as the solution to a Poisson 
equation, driven by the continuity equation 13, in the form 


l($) 



(24) 


Finite Element Penalty Algorithm 

As a consequence of the 3DPNS equation set completion, equations 15, 
16+P3, 1 7 -i 19, 21 and 24, define a well-posed, i ni tia 1 -boundary value 

statement for the dependent variable set q^(x^) -*• (q(x j .)} = {puj,pu £ ,pH, 
pk, p E , p c , p p , *1. The equation of state, p = p(p,H), and equation 20 are 
algebraic definitions for the remaining seven members {p, ulu'. }. 
Therefore, the first nine members of { q ( x ^ ) } are eligible for constraint 
on the solution domain boundary, 3Q = 3R x Xj, by a linear combination of 
Dirichlet and Neumann boundary conditions of the form 


*<V E 


et 
^ l 




a 2 


d \ 

3x„ 


A 01 

a 3 


= 0 


(25) 
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In equation 25, the al are defined to enforce the appropriate constraint 
for each variable, c.f.. Baker [2, 12]. Since the remaining seven members 
of (q(x.)} are defined by algebraic equations, no boundary conditions are 
appropriate. Finally, the first six members of (q ( x ^ ) } are required defined 
on the initial solution plane, = R 2 x xj(0), by an appropriate initial 
condition, {q Q (x £ , Xj(0))}. 


The complete derivation of the finite element penalty constraint 
numerical solution algorithm is given in [1, 12], Briefly, the semi- 
discrete approximation for each member of the set q^x^.) is formed by the 

union of elemental approximations q e (x.) as, 

J 




‘i'V 




(26) 


In turn, each elemental semi-discrete approximation, valid on the 
representative finite element domain = R« x Xi, is formed as an expansion 
on the cardinal basis {N^(x £ )}, the members of which are (typically) 
polynomials complete to degree k, in the form 


q a (x j ) E { W }T { Q ! ( x i)> e (27) 

In equation 27, {•} denotes a column matrix, superscript T its transpose, 
subscript e denotes pertaining to R 2 , and the elements of f QI } g are the 
evaluation of the semi-discrete approximation at the nodal coordinates of R|. 

A basic requirement in any algorithm construction is a formal statement 
regarding constraint on the error formed by employing the semi-discrete 
approximation for the differential equation set. The finite element 
algorithm construction requires the semi-discrete approximation error to 

be orthogonal to the basis employed to construct q^(x.). For all members 

h h J 

of {q } except u £ , the resultant error constraint statement is. 


(N k (x ? )} L( q^)dx + 3 (N k (x £ ) H(q£)dx = {0} 
R 2 ' 3R 


(28) 
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where 6 is a scalar multiplier selected to achieve cancellation of the 
middle term in equation 25. The error extremization statement for the 

h U 

members u" of {q n } is. 


tN kW> } 


L <“J> 




«!>] 


dx +e|{N k (x A )H(uJ)dx 


3R 


+ X 

R 2 


3{N k } 


3x, 


L(^ h )dx = { 0 } 


(29) 


where A is an arbitrary parameter modifying the penalty term which 
constrains the error extremization by the continuity equation (error). 

Equations 28-29 define the finite element penalty algorithm semi-discrete 
error constraint statement for the 3DPNS equation set. For the non-initial - 
valued dependent variables, equation 28 yields the linear algebra statement 
(FI) = {0} , recall equation 1. For the initial-valued variables, equations 
28-29 yield a coupled ordinary differential equation set, 

[A] + {B(.QI ) } = {0} (39) 

which is transformed to a linear algebra statement using a Taylor series, 
for example, 

(FI) = (QI} j+1 - (QUj - Ax, fQI }j +() + ••• - {0} (31) 

where superscript prime denotes the ordinary derivative and 0 > 0 implies 
an implicit statement since equation 30 is quite nonlinear. 

Hence, the final fully discrete approximation error constraint statement 
is the nonlinear algebraic equation set, 

{FI (k. A, 6, Axl {QI }) } = {0} (32) 

where 1 < I _< 16, see equation 1. The Newton algorithm solution for 
equation 32 is given in equations 2-4. 
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DISCUSSION AND RESULTS 


Generalized Coordinates 

A basic requirement for a 3DPNS algorithm is geometric versatility, such 
that the discretization of the transverse plane solution domain R 2 can be 
efficiently embedded within a boundary comprised of the union of select 
aerodynamic surfaces and freestream interfaces. The term "generalized 
coordinates" has gained acceptance in describing an algorithm construction 
suitable for use with a regularizing, boundary fitted coordinate transfor- 
mation. This construction also impacts directly on the efficiency of an 
implicit algorithm, since it usually facilitates factorization of the 
linear algebra statement Jacobian, recall equation 2. An objective of the 
current project is to construct and evaluate a matrix tensor product 
approximation to the Newton Jacobian of the finite element 3DPNS algorithm [8]. 

The required step is to derive the generalized coordinates form of the 
finite element penalty algorithm, equations 28-29. A multitude of procedures 
are available to generate regularizing transformations, c.f., [13], and each 
may be viewed as generating an approximation to the mapping, 

x i = x i ( n j ) (33) 

at a finite number of coordinate triples on the domain ft = R 2 x Xj. For the 

space-marched 3DPNS equation set, equation 33 may be conveniently decomposed 

into an Xj-oriented grid-stretching transformation, X]= x^.), and a 

J 

regularizing transformation on R 2 mapping the boundaries 9R onto coordinate 
surfaces of j = 1, 2. Denote the generated set of coordinate pairs on 
R 2 as {XI} , I = 2, 3. The number of entries in {XI} equals twice the mesh 
characterization of the discretization UR| of R 2 . (For example, for a 41X41 
mesh, there are 2(41)(41) = 3362 entries in {XI}.) On any (each) finite 
element domain R 2 , the specific form of equation 33 is, 

x i = (N k (n) } T {XI } e (34) 

for x. c R 2 . In equation 34, the elements of { XI } g are the appropriate 
members of {XI}, and the interpolation basis {N^(n)} is (potentially) 
identical to the approximation subspace for q!?, equation 27. 
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Equation 34 is of general utility since the elements of {N^n) are well 
known for kj> 1, using either triangular or quadrilateral shaped finite 
elements R*, with or without curved sides. The algorithm requirement is to 
transform the derivatives 9/9x! and 3/3x A , see equations 15-19, 21-24, as 
they appear in the penalty algorithm statement, equations 28-29. A grid 
stretching parallel to the Xi coordinate direction introduces additional 
derivatives. i.n the x^ plane, upon transformation to the t- coordinate system 
of the form [8]. 


3xj 3 1 



- h 3 


_ 2 _ = 3 _ 
9x 3 a t 



(35) 


Therefore, the general form of the ini tial -val ued partial differential 
equations 15, 16+23, 17-19, of the 3DPNS set is 


L <«a> * & 


u ,q a j + 


3x ( 




+ s =0 


ci 


(36) 


Table 1 lists q^, t and s a for 1 £ a _< 6 in equation 36, and 2 £ l <_ 3. 


Table 1. Variables and Parameters in Equation 36 
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The generalized coordinates construction requires transformation of the 
planar divergence operator 3/3x^, when equation 36 is inserted into equation 
28. Using a Green-Gauss form of the divergence theorem on the second term 
in equation 36, yields 


M k } L 


'a 


dx + 


R 2 
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{N k H(q£)d0 = 
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(N> 

'3 

TT 

h h) 

rvl 
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dx 


+ A (N) 


3R 
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n £ d0 


3{N) 9n k 
3n k 3x £ 


f, . 1 

h h . h 

l 1 «J 
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dx 
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ah . a 3q a 9n k p + a « 

q a + 32 InjT" n £ a3 . 


da 


3R 


(37) 


The evaluation of equation 37 is accomplished by performing the integrals 
on an element-by-element basis, and assembling the resultant contributions 
into {FI} using the matrix assembly operator S e » c.f., [12, Ch. 2]. Since 
the elements of {N^(n)} are known functions of n, the only requirement is 
evaluation of 3ri( ( ,/3x^ on a generic finite element domain R 2 . From equation 
34, the 2X2 square matrix defining the Jacobian of the forward transformation 
is 
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(38) 


[j] e = CJ({xn e )] =- 



Thus, the elements of the inverse transformation Jacobian are. 




1 


det J, 


[Cl 


(39) 


where [C] is the 2X2 transformed cofactor matrix of [J] g , the elements of 
which are algebraic functions of n k and {XI } £ . The differential element dx 
in equation 37 becomes 

dx = det J g dn (40) 


Finally, it is convenient to define the contravariant components of the 
convection velocity semi-discrete approximation u^ on as 


u e k = det J e 





(l-h £ ) CC] k£ u £ 


(41) 


With equations 38-41, and recasting evaluation of equation 37 as the 
assembly of integrals over UR| , the generalized coordinates form of the 
finite element algorithm statement, equation 28, becomes 



det J e {N} {N} 1 


{U1 QI } e + (SI} e dn 


R 


2 
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{UBARKlI {N} ~ {N} fN} 1 
e df)^ 


{QD e dn 


R 2 


{ETAKL}g {N} {N} {N} T {TAUIL) e dn 

k 


'I 


+ ic 


{N} [a i q 4 + 23 3 det J e do 


9R 


+ 0 


{N} [(l-h £ ) OIL} 1 {N} {N} T {QI} ( 


3R 


+ UO T {TAUIL} e n L det J c jdo 


(42) 


In writing equation 42, it has been assumed the only variable possessing 
a constrained normal derivative boundary condition is stagnation enthalpy. 
Hence, the fourth term contains only q^, and the normal derivative term in 
equation 25 is cancelled by the cor res pons ing term in the third integral in 
equation 38 by defining B = k, see Table 1. A second assumption is that 
det is adequately represented as an elemental scalar, which is a commission 
of interpolation error (only) on a sufficiently refined mesh. The elements 
of the direction cosine matrix {ETAKI > e are defined by the interpolation. 

— — {N. } T {ETAKL} (43) 

det J e 


3n ? 

l*k 
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For example, defining k = 1 in equation 43, for the bilinear tensor 
product basis on straight-sided quadrilaterals, numbering nodes counter 
clockwise, and defining the elements of the nodal array {XI > e as {YJ,.ZJ, 
1 1 J 1 ’ the four arra y s 2( ETAKL) e are, c.f., [12, Ch. 8]. 


2 { ETAKL } g 


( " 
Z4-Z1 

- 

•> 

Y1-Y4 

12-12 


Y2-Y3 

12-12 

< 

* < 

Y2-Y3 

Z4-Z1 


Y1-Y4 

. (2,2) , 


. ( 2 . 3 ) . 


j \ j 

e e 


Z1-Z2 


Y2-Y1 

Z1-Z2 


Y2-Y1 

Z4-Z3 

. 

Y3-Y4 

Z4-Z3 


Y3-Y4 

(3,2) 

1 

l (3,3) > 


e e 


(44) 


The indices in the parenthesis in equation 44 denote (K, L). Finally, the 
matrix elements of { UB ARK } e are defined by the interpolation. 


U k = {N k } T (UBARK) e 


(45) 


Proceeding through the algebra yields 

{UBARK) e = {ETAKL) e {N k } T {UL> e (46) 

Since the evaluation of equation 46 is at nodes, the elements of {N. } reduce 

T * 

to the Kronecker delta; hence, {N k } { UL } e simply selects the correct nodal 
value out of the elemental arrays {UL} g , L = 2, 3. 

In equations 42-46, the indices K, L are tensor summation indices, while I 

L 

denotes the appropriate member of {q k } . Every matrix denoted by a subscript 
e in equation 42 is independent of the n coordinate system spanning R*, hence 
can be extracted outside the integral. The expressions remaining within the 
integrand are polynomials on q , for all definitions of completeness k of the 
semi-discrete approximation, recall equations 26-27. The order of {N k } depends 
upon k, as do all the integrels in equation 42, and each is readily evaluated 
using numerical quadrature. The basic structure of equation 42 is thus defined, 
using a standardized hypermatrix nomenclature [12], as 
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= S e det J e [B200] -^-{UIQI } g + det J 0 [B200] {SI > e 

- {UBARKig [B3QK0] - 1-h^ n K det J g {UK}J [A3000] {QI } g 

- {ETAKL}T [B30K0] - n L det J g [A200] {TAUI L } e 

+ k 6 i4 det J fa 4 [A200] {QI} + a 4 {A10}] e {0} (47) 

Note that equation 47 is in the form of equation 30, with the definition 
[A] = S g [det J e [B200]]. The generalized coordinates finite element penalty 
algorithm involves definition of only three distinct standard matrices, 

[B200] and [B30K0], K = 2, 3, on R 2 These are independent of R 2 and depend 
only on the degree k of the semi-discrete approximation. In addition, there 
are the three standard matrices [A200], [A3000] and {A10}, defined on the 
boundary of R 2 , hence evaluated on a one-dimensional element R*. In the 
second and third lines of equation 47, the terms involving [A ••] are the 
residuals from use of the Green-Gauss theorem, see equation 37. The terms in 
the fourth line result from the heat convection boundary condition for 
stagnation enthalpy. The tensor indices range 2 <_ (K , L) <_ 3, the subscript 
J< takes on the value of K , and n^ are unit normal vectors to R 2 , and 
1 ^ I _< 16 denotes the ordered members of the dependent variable array. 
Finally, det J g is the average value on R 2 of the determinant of the transfor- 
mation Jacobian, equation 38, and { ETAKL } e and { UBARK } e are element dependent 
matrices defined in equations 43-46. The Appendix lists the defined matrices 
[A ♦ •] and [B* •]• f° r k = 1 in equation 27. 
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Tensor Matrix Product Jacobian 

The form of equation 47 is well suited to formation of a tensor (outer) 
matrix product approximation to the Newton iteration algorithm Jacobian, 
see equations 2 and 4, upon definition and use of a tensor product cardinal 
basis {N^}, equation 34, on quadrilateral finite element domains (which 
could be curved-sided , k > 1). in this instance, the Newton algorithm 
Jacobian is approximated by the tensor matrix product construction, 

[J( {QI })]=> [J 2 ] 0 [J 3 ] (48) 

where 0 denotes the outer product [14]. Equation 2 then becomes of the form 

[J 2 3 ® [J 3 ] (SQU = - (FI) (49) 

Making the definition [J 3 ] { 6QI } = { 6 P I } , equation 49 is solved by sweeping 

the mesh parallel to n 2 » to determine { 5PI } , and then sweeping the mesh 
parallel to n 3 to complete the solution, i.e., 

[J 2 ] UPI) e = - {F I } e 

M (5QI} e = (6PI } e (50) 

The attraction of the approximate solution statement, equation 50, is 
replacement of the large sparse matrix [J], equation 4, by two block 
tri- or penta-diagonal matrices [J ], with the corresponding significant 
reduction in computer storage and LU decomposition computer CPU time. 

The detraction of the statement is degradation of the convergence rate 
associated with the exact Newton algorithm statement. Of course, for the 
16-dependent variable 3DPNS equation set it is impossible to construct and 
use the exact Jacobian. Hence, the trade-off is degree of approximation 
versus associated cost. The critical measure is convergence rate since 
this determines the number of passes through the linear algebra statement. 

Study Problem Statement 

The tensor matrix product approximation to the 3DPNS Newton Algorithm 
Jacobian is formed by evaluation of the appropriate integrals in !){FI}/3{QJ} 
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on a one-dimensional domain R*. The CMC:3DPNS code [10] was reorganized to 
permit construction of a one-dimensional Jacobian, independent of the sub- 
routines used to generate {FI}. The basic issue is convergence, and the 
study problem is solution of the 2DPNS equations for laminar and/or turbulent 
flow. The code was indexed to permit use of the original algorithm formula- 
tion, as well as several variants of the tensor product construction. The 
2DPNS algorithm construction for turbulent flow using the k-e closure with 
a Reynolds stress algebraic model was established, however, the detailed 
discussion is limited to the equation system for laminar isoensrgetic flow. 


The governing 2DPNS differential equation system, see equation 36, for 
{ q> = {u, v, p, <(>}, is 


L(u ) = |j- (uu) + p 1 



V 

Re 



(51) 


L 6 (v) 




vv + 


£_ . v_ 

p Re 


DV_ 

ay 


= o 


(52) 


L(p) 


a_ i_ aji 

ay p ay 



+ v 


Dv 

ay 


1_ D_ f 3v' 
Re 9y ( v 9y J 


= 0 


(53) 


LU) 

ay 2 


Du Dv 
Dx ” Dy 


(54) 


In equation 51, p' = ^ ^ is an input parameter (corresponding to the 
p r (x l5 x ) solution, equations 21-22, for 3DPNS). Inserting equations 51-54 
into the finite element penalty algorithm, equations 28-29, see equation 42, 
proceeding through the algebra, equation 47, and inserting the results for 
u n and v n into the trapezoidal rule (S-j) form of equation 31, equation 32 
becomes the set. 


{FI} = S { FU } 
e e 


S 

e 




U.}' [A3000] 
J ® 





-y- {V)J [A3001 ] {U) e + A e {P') e + {NU} T [A301 1 ] { U } 

e 


(55) 


J + l ,j 
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{F2} E S {FVI = S 


? VU $ +1 + Uj}^ [A3000] 


{V }P 

lV 'j+l 


(V h 


Ax 


(V}J [A3001 ]{V) e + 1 [A201]{P) e 


+ 7 - {NU}p [A301 1 ] {V } 
e e 


j+1 , j 


+ |TUJ e [A 201 ] {Sd>} £ 


j+1 


(56) 


{ F3 } E 


S {FPL 
e e 


^-[A211]{P} e + 


m] [A3 01 0] {VP} ( 


+ J— {V}J [A3011] {V} g 


EuA 


J j+1 


(57) 


{ F4 } 


S{F<J>L = S' 


-[A21 1 ] {cj>} + A [A200] { UP } + [A201 ] {V}, 


(58) 


’j+1 


The newly defined standard matrices [A* •] are also listed in the Appendix. 
Further, det J g e A g , the measure of R^, and its occurrence has been cleared 
throughout equations 55-58, such that [A***]are integer arrays with a common 
divisor. The continuity equation has been explicitly employed to clear extra 
terms in equations 55-56, {«}j + -j denotes evaluation with the current iterate 
at x i+i> (•}. denotes the converged solution at the previous step (x.), 

indicates evaluation with both { • } j + i and { # }j, followed by addition, 
and [*]j + i denotes evaluation with { * > j + i • The arrays {UP > e and {VP} g contain 
a finite difference approximation for ^ {U) e and {V> e , evaluated at 
x j + 1 • The elements of { NU} e are v/Re, where Re is the Reynolds number, Eu is 
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the Euler number, and {P'l contains — 4^. 

e p dx 

The last term in equation 56 is the continuity constraint penalty term, 

the selected form of the X- modified integral in equation 29. is a constant 

of order unity, and fUJ e is a diagonal matrix with elements taken from {U} g 

in the same order. (This multiplier improves the resolution of {V} in the 

lower reaches of the boundary layer. The penalty term dependent variable is 
{S4>) e , defined as the sum of the previous solutions to equation 58 at step 
x j + i » 1 • G - 5 


{Sep} 


P+1 

j+1 


MS*>P + , ♦ 

P+1 

{4>} j+l 


P 

f 

p 

P+1' 

= (S<J>}j + -| + 


+ 

> 

P-1 

P-1 

P 

= (S<J>} j+1 + 

W J+ i 

+ {64>}j + -] 


+ 



P+1 

+ W j+1 


, etc. 


(59) 


Furthermore, {S4> ) j + i = which is the identification as well for the 

first estimate of the dependent variable set, i.e. f {QI {j + i = This 

definition of (S<f>}, as the penalty variable, is a considerable modification 
of the orginal algorithm construction in CMC:3DPNS, to permit solution for 
{ <5<p } within equation 2, rather than a separate Poisson solution for {<{>}. 


Equations 55-59 constitute the linear algebra specification of the 2DPNS 
algorithm. The Newton algorithm Jacobian matrices are derived from the 
definition, equation 4, using equations 55-58. They are identical to the 
factor [J^] of the tensor product approximation for the 3DPNS statement. 
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equation 50. Recalling that [J] = S e C J J e , the elemental components of the 
4-block, tridiagonal Jacobian .{J 2 ] are: 

[JUU] e = A e {U}g [A3000] + T{V}g [A3001 ] + J- {NU}^ [A3011]] 

[JUV] e = [A3100] 

[JUP] e = [0] 

iml e = CO] (60) 


[JVV] e = ±A e {U 


P 

j+1 


+ u. 

J 



l [A3000] 

(V}J [ [A3001 ] + [A3100] 


j- {NU}g [A3011] 


[JVU] e = ^A e (Vj +1 - V.}] [A3000] 

A e 

CJVP] e = [A201 ] 

[JV<D] e = rUJ e [A201] ( 61 ) 


[JPP] e [A2H] 

H e 

[ JPU] e = J- {VP)J [A3010] 

t JPV 3 e » -n {U>e [A301 0] + (V}T 

e 

[JP<D] e = CO] 


f \ 

[A3011] + [A3110] 

. 


(62) 
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[J#] e = 4- [A211] 
e 

[J<j>U] e = aA e [A200] 


[J4>V] e = [A201 ] 


[J4>P] e = 


[ 0 ] 


( 63 ) 


In equations 60-63 all evaluations are made using { - } j + -j unless otherwise 


noted. The parameter "a" in equations 62-63 denotes the fraction of (QI 
used in the finite difference approximation for ^.{QI} at x^ . 

Numerical Results 


j+1 


Numerical experiments have been executed to evaluate the 2DPNS algorithm 
using various approximations to the Newton algorithm Jacobian, equations 
60-63. The exact, construction is 

P 


S 


e 


[JUU] e 

[JUV] e 

0 

0 

[JVU] e 

[JW] e 

[JVP] e 

[jv*] e 

[JPU] e 

[JPV] e 

[JPP] e 

0 

[J<t>U] e 

[J4>V] e 

0 

[ml e 


= - fFI}P +1 (64) 


j+1 


For reference, both the 2D and 3D algorithm constructions currently oper- 
ational in the CMC:3DPNS code approximate equation 64 by its diagonal 
entries only. This decouples the dependent variable set solution in the 
nondiagonal terms premultiplied by Ax, eauations 60-63, but retains the 
Ax-coupling through convection, pressure gradient and viscous-turbulent 
effects. Further, the original algorithm Poisson equation solutions are on 
{p} and {<)>}, rather than {6P} and {6<j>}, and both are computed using the most 
recent solutions {U}?^, {V}?*}, etc., rather than (QI )j +1 , see equation 2. 


The 2D algorithm construction for solution of {P} and {(f>) was rearranged 
and the code modified to permit PNS solutions using various approximations 
to [J] . The nearest equivalent to the original algorithm construction is. 
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p 


[JUU] 

[JW] 

[JPP] 

where {FP} and {Ft};} are now evaluated using {QI }P +1 . Defining the penalty 
variable as {S<J>}, see equations 56 and 59, reproduced the original algorithm 
solution data to within negligible differences. The (finally) selected 
laminar flow standard test case used a geometrically nonuniform, unstretched 
grid of M=56 elements spanning 36 q . The initialization of {U} Q employed the 
Blasius solution, and the boundary layer thickness doubled over the integra- 
tion length Ax = 4, achieved in nominally 20 nonuniform integration steps Ax 
at Re = 3.2x1 0^. 

Table 2 summarizes data confirming the modified 2DPNS algorithm can 
generate an acceptably accurate solution on the M=56 mesh, with nodal coor- 
dinates as listed in column 1. Columns 2-4 compare the axial velocity 
profiles of the Blasius solution, the direct boundary layer (2DBL) solution 
using [JUU] and a direct trapezoidal rule integration of the continuity 
equation 14, and the basic 2DPNS algorithm. The matrix iteration convergence 
requirement was |6U| < 10 5 , and peculiarities in starting the 2DPNS solution 
were exactly mimicked in generating the 2DBL solution. Hence, both numerical 
solutions are displaced axially from the Blasius solution by a negligible 
increment. The agreement between the 2DBL and 2DPNS solutions is excellent in 
both u 2 and in u 2 , see columns 5-6 of Table 2. The transverse velocity 
distribution is the sensitive measure, since it ranges over five digits of 
significance over 6 (see arrows) and is required to implicitly approximate a 
vanishing normal derivative at y/S = 0. 

With this confirmation, the principal interest is efficiency. Table 3 
summarizes comparison data on convergence of the approximate Newton algorithm 
Jacobian constructions. All computations are done in single precision using 
a 32-bit word. Using the Blasius solution for (V(x,y)} the 2DBL data of the 
second column indicates use of [JUU] yields approximate quadratic convergence 
to E(-7) for a scalar equation. Appending the direct continuity equation solution 


{6QIK 


P+1 _ 


j+1 




(65) 




Table 2 


Accuracy Comparisons, Laminar Boundary Layer Test Case, AX = 4. 


NODES 

J j 

AXIAL 

VELOCITY - u-j 

X 10" 1 

TRANSVERSE VELOCITY 

- u 2 X 10 3 

2y/106 o 

Blasius 

2DBL 

2DPNS 


2DBL 


2DPNS 

0759$ 

099568 

099965 

099966 


215451 


216400 

07777 

099950 

055945 

055546 


215312 


2 l 6260 

07555 

059523 

099916 

0559 17 


215106 


216054 

07333 

099883 

099872 

099874 


214807 


215755 

07111 

059624 

055807 

055810 


214379 


215325 

06688 

059 74 L 

099714 

055720 


2 13778 


214731 

06666 

099622 

099584 

055551 


212949 


213908 

06 44 4 

059453 

055401 

099413 


211827 


212754 

06222 

099221 

099153 

099169 


210333 


211314 

05559 

098508 

09 6818 

098840 


206362 


209379 

057 7 7 

058454 

098375 

098405 


205879 


206657 

05555 

097555 

097798 

097838 


202725 


203771 

05333 

097264 

097060 

097111 


198825 


199903 

05111 

056355 

096130 

C96 195 


194067 


195159 

04666 

095314 

094977 

055057 


188438 


185588 

04666 

053594 

093570 

093666 


181822 


183013 

04444 

052390 

09 1 8 60 

09 1994 


1742 15 


175432 

04222 

090458 

089881 

0900 12 


165623 


166868 

03595 

088196 

087551 

087699 


156092 


157361 

03 77 7 

065588 ' 

0648 74 

085037 


1457C4 


146962 

03555 

08262 1 

03 1841 

0820 19 


134579 


135835 

03233 

079293 

0 7 84 52 

078640 


122871 


124117 

03111 

075605 

074713 

0749C8 


1 10761 



111928 

02 68 6 

071565 

070637 

070836 


058450 


059573 

02666 

067200 

066244 

066443 


086147 


087183 

02444 

062524 

061563 

061758 


074066 


074966 

02222 

057565 

0566 22 

056808 


062411 


063274 

01595 

052339 

051455 

051630 


051374 


051536 

01722 

045546 

044732 

044887 


038712 


039737 

01480 

039421 

038696 

038832 


0289 13 


025649 

01270 

033 568 

03334U 

033457 


02 1443 


022085 

0 l C 6 7 

025 160 

C28617 

028 7 1 7 


015761 


016265 

00928 

024551 

024483 

024568 


011551 

^ 

011663 

00790 

021266 

020858 

020931 


008393 


0007 15 

C06 7C 

0 18C52 

017704 

017766 


006037 


006337 

00566 

0 15245 

014953 

015005 


004316 


0045 10 

i 00475 

012803 

012556 

012599 


003045 


003183 

i' 00396 

0 106 76 

010472 

C 10506 


002116 


002250 

00327 

008627 

008661 

008691 


00 1458 

^ 

001577 

00267 

007219 

007080 

007105 


000988 


001074 

00216 

005620 

; 005708 

005728 


000649 


000701 

00170 

004604 

i 0C4514 

004529 


000409 


000436 

00 13 1 

003546 

' 0034 76 

0034 6 8 


000243 


000258 

00C57 

002626 

0025 74 

002583 


000134 



000146 

00067 

001626 

001750 

001796 


C0C065 


000079 

00C4 1 

001130 

001108 

001112 


000025 


000040 

OOCIS 

000525 

000515 

000517 


000005 


0000 16 

0 

0 

— 

0 

0 


C 


0 
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Table 3 


Newton Algorithm Convergence Summary 



Maximum Residual 

Iteration 
Index (p) ( 

2DBL Solution 

2DPNS 

Solution 

Blasius {V} 

Computed {V} 

Original 

Revised 

1 

-1.2(-2) 

-1.2(-2) 

-1.2(-2) 

- 1 • 2 ( -2 ) 

2 

1 . 4 ( - 4 ) 

5.4(-5) 

1 . 4 ( -4 ) 

1.5(-5) 

3 

5.8(-7) 

1 . 3 ( -5 ) 

-7. 6(-5) 

-9.4(-5) 

4 


3 . 5 ( - 6 ) 

4.2(-5) 

-5.3(-5) 

5 



-2 . 3( -5) 

1 . 8 ( - 5 ) 

6 



l.K-5) 

4 . 7 ( -5 ) 

7 



8.4(-6) 

5. 1 (-5) 

8 




3.2C-5) 

9 




5.7(-6) 


for the 2DBL column 3, the convergence is quadratic only to E ( -4 ) for 

the two equation system. This decline in convergence rate for | 6U I < 10"^ 
characterizes all the 2DPNS algorithm solution constructions as well. The 
fourth column in Table 3 summarizes convergence of the original PNS algorithm 
diagonal [JQQ], and solution for {P}?*] and using {QI}?*], which is 

monotonic and approximately linear for 1 6U j < 10 - 4. The last column gives 
these data for the new algorithm construction, equation 64, with [JPU] and 
[J4>U] omitted (since they cause instability and eventual divergence). 
Convergence is quadratic to E(-4), and thereafter is nonmonotonic and subl inear 
(although these are extremum residuals occurring at different nodal coordinates 
at any iteration) . 
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The current practice in CMC:3DPNS is 4-5 iterations per step with conver- 
gence set at l<$Q| max < 10 - ^. Over this range the convergence character of 
the old diagonal [JQQ] construction and the coupled construction is nominally 
identical. Additional examinations were conducted to assess reasons for the 
poor convergence rate below E(-4). Since the action of the penalty term is 
applied modulo a discrete (second-order) approximation to 9/3y, the {V} 
solutions generated at ten iterations/step or more will eventually exhibit 
"wiggles." The current practice is to element average either {6V} or {V} 
when this occurs. In Table 4, column 2 summarizes the standard test solution 
for {V} using the original algorithm (for reference), column 3 contains the 
revised algorithm data using (V) averaging, and column 4 contains the same 
algorithm solution without {V} or {<5V} averaging. Close examination of the 
data in column 4 verifies periodic occurrences of local flat spots which will 
eventually grow into a 2Ay wave. Since a {V} and/or (6V) average is tantamount 
to not using the coupled implicit solution vector {SQI} as computed, this 
operation could contribute to poor convergence. A numerical test verified this 
to a limited extent* compare column 2 of Table 5 to column 5 of Table 3. 

Alternatively, analysis determined that the "wiggles" in (V) can be traced 
back directly to {<f>}, hence {S<t>}, which if smoothed prior to use in equation 
56 would not generate the discrete wave. The correct way to obliterate a 2Ay 
wave is to use a Shuman-type digital filter, cf., [12, Ch. 4]. As an approx- 
imation, {S4>}j^| was simply averaged. The accuracy of the resultant solution 
for {V}, column 5 in Table 4, is nominally identical to the (V)-averaged 
solution data, except directly adjacent to the wall where the averaging removed 
the critical sensitivity. This operation did improve the iteration convergence, 
column 3 in Table 5, mostly in returning it towards monotonicity below 

| <5U I < 10' 4 . 

1 'max 

The Newton Jacobian for these tests remained incomplete, since including 
either [JPU] or [J<f>U] would destabilize the algorithm. (Both these terms involve 
a discrete approximation to an axial derivative {QIK, see equations 62-63, 
to which the PNS penalty algorithm is quite sensitive.) The semi-implicit 
evaluation of {QIK, of the form. 
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Table 4 


2DPNS Algorithm Accuracy 


NODES 

2y/l06 

0 

JUU.JVV 
PPRESS 
Q‘ Imp. 

J(I,J) 
V Avg. 
Q’ Imp. 

07999 

216400 

214784 

07777 

216260 

214629 

07555 

216054 

214402 

07333 

j 215755 

214076 

07111 

215325 

213613 

0668 8 

! 214731 

212970 

06666 

213908 

212090 

06444 

212794 

210907 

06222 

211314 

209346 

05999 

209379 

207320 

05777 

206857 

204738 

05555 

203771 

201506 

' 05333 

199903 

197531 

05111 

195159 

192727 

0486 6 

189588 

187026 

04666 

183013 

180377 

04444 

175432 

172760 

04222 

166868 

164184 

03999 

157361 

154695 

03777 

146962 

144378 

03555 

135835 

133351 

03333 

124117 

121764 

03111 

111928 

IG5754 

02688 

099573 

097631 

02666. 

087183 

085485 

02444 

074966 

073566 

02222 

063274 

061926 

01595 

051536 

050403 

01722 

039737 

039301 

01 480 

025649 

029607 

01270 

022085 

021952 

01087 

016265 

016171 

00928 

011883 

01 1863 

00790 

0087 15 

008671 

00670 

006337 

0063C2 

00566 

0045 10 

004551 

00475 

003183 

003253 

00356 

002250 

: 002288 

00327 

CO 15 77 

001575 

C 026 7 

001074 i 

001058 

002 16 

000701 

000655 

CO 170 

000436 

000448 

00131 

000258 

0002 es 

00057 

000146 

000180 

0006 7 

000075 

000112 

0004 1 

000040 

000064 

0CC19 

000016 

000029 

0 

- 

0 

0 


Comparisons 

, Laminar Test 

Case 

J(I,J) 

J(I,J) 

J(I,J) 

No V Avg. 

S<J> Avg. 

S<j> Avg 

Q’ Imp. 

Q 1 Imp. Q 1 

Sem-Imp 

220248 

218018 

220856 

220080 

217866 

220733 

219993 

217643 

220549 

219559 

217321 

220280 

219416 

216865 

219893 

218562 

216229 

219347 

218138 

215356 

218586 

216796 

214180 

217550 

; 215631 

212625 

216159 

; 213703 

210605 

214330 

i 211318 

208022 

211964 

208441 

204779 

208961 

! 204574 

200785 

205219 

| 200052 

195950 

200644 

194729 

1902 13 

195140 

187884 

183443 

188642 

180835 

175760 

181157 

172174 

167138 

172609 

162172 

157274 

162973 

152939 

147017 

152646 

140295 

135736 

141402 

129053 

123684 

128928 

117556 

111736 

117288 

102326 

055510 

104491 

092874 

086025 

089834 

078062 

C757C9 

081177 

065551 

06 1870 

064422 

057151 

052053 

055122 

038727 

035070 

042943 

033067 

026164 

028567 

022816 

022300 

024686 

014982 

015474 

016477 

014978 

011332 

011714 

008086 

008553 

009200 

004866 

006070 

006634 

006336 

004428 

004565 

003469 

003071 

003123 

000324 

002187 

002038 

001338 

001460 

001242 

002293 

001056 

000790 

001278 

0C06C5 

000633 

000155 

000582 

000575 

000164 

000484 

000406 

000651 

CCC038 

000072 

000843 

-000280 

000305 

000503 

-000419 

A 

000526 

000058 

-000325 -000435 

0 

0 

n 
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TABLE 5 


Newton Algorithm Convergence Summary 
Revised 2DPNS Algorithm 



Maximum Residual 



Semi-Impl icit (QK 

Iteration 
Index (p) 

Implicit {Q}' 
No {V} Avg. 

Implicit {Q}' 
{Sd>} Avg. 

{S4>} Avg. 
Full [J] 

(Sd> } Avg. 
Orig. [J] 

1 

-1.2(-2) 

-1.2(-2) 

“1 - 2 ( -2 ) 

-1.2(-2) 

2 

-8 . 0(-5) 

2 - 2 ( - 5 ) 

-8.2(-5) 

9.4(-5) 

3 

-5.0(-5) 

3 . 2 ( - 5 ) 

- 3 . 3 ( - 5 ) 

-4.0(-5) 

4 

3. 6(-5) 

3 . 4 ( - 5 ) 

5.0(-5) 

-2.8(-6) 

5 

4.K-5) 

2.2(-5) 

4.1(-5) 


6 

2 . 6( - 5 ) 

7 . 5 ( -6 ) 

2 . 0( - 5 ) 


7 

1 - 4 ( - 5 ) 


8 . 7 ( -6 ) 


8 

-l.l(-5) 




9 

-1 . 3 ( - 5 ) 


i 






( 66 ) 


Is the standard procedure in CMC:3DPNS. With this simplification, the 
parameter "a" in equations 62 and 63 is zero. As a consequence, [J<pU] 
vanishes identically and [JPU] can be included without destabilization. 
Using {S4>} averaging, the accuracy of this algorithm form is nominally 
unchanged, see column 6 in Table 4. The resultant Newton Jacobian is 
exact and these results are a very modest improvement in monotonicity 
of convergence, see column 4 of Table 5. This convergence character is 
the closest to the original algorithm, recall column 4, Table 3. Of more 
significance, insertion of {S<J>} averaging into the original algorithm 
construction significantly improves convergence, see column 5 of Table 5. 
This operation yields the PNS algorithm as efficient as the direct 2DBL 
solution, recall column 3 of Table 3. 



SUMMARY AND CONCLUSIONS 


A generalized coordinates form of the penalty finite element algorithm for 
the 3-dimensional parabolic Navier-Stokes equations for turbulent subsonic flows 
has been derived. This algorithm formulation requires only three distinct hyper- 
matrices in its formulation and is applicable using any boundary fitted coordi- 
nate transformation procedure. The tensor matrix product approximation to the 
Jacobian of the Newton linear algebra matrix statement has been derived. The 
Newton algorithm has been restructured to replace the large sparse matrix solution 
procedure with a grid sweeping procedure using a-block tridiagonal matrices where 
a equals the number of dependent variables. 

The principal purpose of the reformulation is to improve solution economy. 

With the restructured Jacobian, solution economy is linearly dependent on the 
convergence (rate) of the Newton algorithm. A series of numerical experiments 
were performed to evaluate convergence as a function of Jacobian completeness and 
off-diagonal coupling. The results of these studies indicate that the favorable 
Newton quadratic convergence rate is maintained to a residual level of order 10 _ \ 
Thereafter the convergence rate uniformly decreases to linear for residual compu- 
tations in the range 10 _<!> - 10" 5 . Several modifications to the implicitness of 
the algorithm Jacobian and to the overall linear algebra statement were made and 
evaluated. Comparison to exact solutions indicates adequate accuracy is attain- 
able for each of the modifications. 

The results of this study provide required guidance on the appropriate form 
for the tensor product 3DPNS algorithm. The original form of the algorithm, 
employing a diagonal Jacobian, retarded evaluation of the Poisson equation 
solutions, in particular the <f> solution, and {S<J>} averaging for the penalty term 
yields the best Newton convergence performance and solution accuracy for the 
test case. The considerable effort in constructing and coding the off-diagonal 
Jacobian entries appears unrewarded, based upon these data, especially 
considering the added solution costs associated with a block versus scalar tri- 
diagonal matrix. This indication gains considerable importance in the progres- 
sion to turbulent and/or three-dimensional flows, wherein the Reynolds stress 
tensor will almost double the block size. Based on these data, it appears that 
the 3DPNS tensor product algorithm should employ a nominally diagonal tensor 
product Jacobian approximation for the initial-value variables, and should 
retard the Poisson equation solutions, which in themselves can use a scalar 
diagonal tensor product Jacobian approximation. The computation of Reynolds 
stresses would also employ a scalar diagonal form when using an algebraic model. 
This reconstruction of the 3DPNS algorithm would fit directly into the present 
CMC:3DPNS code. It is suggested that this should indeed by the next step. 
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Finite Element Algorithm' Standard Hypermatrices 
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